#!/usr/bin/env bash

####on vérifie que 3 parametres au moins sont passés en argument
if [ $# -lt 2 ];then
	echo -e "\n ERROR : scriptIBD.sh <filename> <mapFile> <OPTIONAL : IBD training sample file> \n" 1>&2
	exit 1
fi

nb=$RANDOM
filename=$1 #récupération du 1er parametre

echo -e "\nDo you want to search ambiguous SNPs in build 132 ? (y/n)"
read choice
if [ "$choice" == "y" ];then
	echo "load snp132"
	# Pour faire le lien avec la build 132
	curl --proxy proxy-upgrade.univ-nantes.prive:3128 -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp132.txt.gz" |gunzip -c |cut -d'	' -f5,18 |awk -F'	' '($2 !=  "1" )' |cut -d'	' -f 1 |sort -k1,1 | uniq |gzip --best > dupsnp.txt.gz 
	echo "see dupsnp.txt.gz"
	echo -e "\nDo you tu remove these SNPs in datas ? (y/n)"
	read choiceRM
	if [ "$choiceRM" == "y" ];then
		gunzip dupsnp.txt.gz
		plink --noweb --bfile $filename --exclude dupsnp.txt --make-bed "new"$dataname
	fi
fi


cut -d" " -f1 $filename.fam  | sort | uniq > ped_$nb # récupération de la liste des familles

dataname=`basename $filename`

#### création des répertoires ped 12 et outplink s'ils n'existent pas déjà
if [ ! -d intermediateFiles ];then
	mkdir intermediateFiles
fi
if [ ! -d intermediateFiles/ped12 ];then
	mkdir intermediateFiles/ped12
fi
if [ ! -d intermediateFiles/outplink ];then
	mkdir intermediateFiles/outplink
fi
if [ ! -d intermediateFiles/forIdcoeff ];then
	mkdir intermediateFiles/forIdcoeff
fi
if [ ! -d intermediateFiles/newtfile ];then
	mkdir intermediateFiles/newtfile
fi
if [ -e $dataname"_new.ped" ];then
	rm $dataname"_new.ped"
fi
if [ -e ibdldreskin_$dataname ];then
	rm ibdldreskin_$dataname
fi

#### pour chaque famille du ".fam" de départ
while read line;do
	
	# récupération des colonnes ped et ind correspondant à la famille 
	if [ "$choiceRM" == "y" ];then
		grep $line $filename.fam | cut -d" " -f1,2 > intermediateFiles/ped12/$line.ped12 
	else	
		grep $line "new"$dataname.fam | cut -d" " -f1,2 > intermediateFiles/ped12/$line.ped12 
	fi

	# création des fichiers d'entrée pour le script1.r
	if [ "$choiceRM" == "y" ];then
		plink --noweb --bfile "new"$dataname --keep intermediateFiles/ped12/$line.ped12 --make-bed --out intermediateFiles/outplink/$line
	else
		plink --noweb --bfile $filename --keep intermediateFiles/ped12/$line.ped12 --make-bed --out intermediateFiles/outplink/$line
	fi

	plink --noweb --bfile intermediateFiles/outplink/$line --transpose --recode --out intermediateFiles/outplink/$line
	
	echo -e "\n########### CREATE FILES FOR idcoeff : $line #############\n"
	Rscript script1.r intermediateFiles/outplink/$line.tfam

	echo -e "\n############# SETTING UP idcoeff : $line ###############\n"
	# Lancement de idcoeff
	./idcoefs -p intermediateFiles/forIdcoeff/pedcoef_$line.p -s intermediateFiles/forIdcoeff/pedcoef_$line.s -o intermediateFiles/reskin_$line

	# modification des fichiers de sortie de idcoeff
	perl modif_idcoefoutput.pl $line intermediateFiles/reskin_$line intermediateFiles/forIdcoeff/corresp_$line.txt

	echo -e "\n########### CREATE FILES FOR IBDLD : $line #############\n"
	Rscript script2.r intermediateFiles/outplink/$line.tped intermediateFiles/outplink/$line.tfam intermediateFiles/forIdcoeff/orderped_$line".p" $line
	perl create_newtfam.pl $line intermediateFiles/forIdcoeff/orderped_$line.p intermediateFiles/outplink/$line.fam
	plink --noweb --tfile intermediateFiles/newtfile/$line"_new" --recode --out $line"_tmpnew"$nb

	# concaténation des données
	cut -d" " -f1-5,7- $line"_tmpnew"$nb".ped" >> $dataname"_new.ped"
	cat ibdldreskin_$line >> ibdldreskin_$dataname

done < ped_$nb

echo -e "\n############# SETTING UP IBDLD ###############\n"
if [ $# == 2 ];then
	./IBDLD -o $dataname -p $dataname"_new.ped" -m $2 -i ibdldreskin_$dataname -method LD-RR -ploci 10  -ibd 3 -hbd  #$2 map file
else
	./IBDLD -o $dataname -p $dataname"_new.ped" -m $2 -i ibdldreskin_$dataname -method LD-RR -ploci 10  -ibd 3 -hbd -unphased $3  #$2 map file, $3  Training Sample File
fi

echo -e "\n############# HANDLE IBDLD OUTPUT (QUALITY CONTROL : kinship(observed) comparison with meanIBD(expected)) ###############\n"
gunzip *.kinship.gz

if [ -f $dataname.kinship ];then
	rm $dataname.kinship
fi

ls $dataname*.kinship > list
touch $dataname.kinship

while read line;do
	cut -d"	" -f1,2,3,4,5,15 $line | awk '{print $1"_"$2"_"$4,$2,$4,$5,$6}' | tr " " "	" >> $dataname.kinship
done < list

cat ibdldreskin_$dataname | awk '{print $1"_"$2"_"$3,$2,$3,$10,$11,$12}' > tmpreskin1_$dataname
cat ibdldreskin_$dataname | awk '{print $1"_"$3"_"$2,$3,$2,$10,$11,$12}' > tmpreskin2_$dataname

Rscript IBDkinship.r $dataname.kinship tmpreskin1_$dataname tmpreskin2_$dataname $dataname

rm list
rm ped_$nb
rm *_tmpnew*

# example de commande ./scriptIBD6.sh datas/brugada7 datas/brugada8.map datas/nvtsfdesir.txt 
# example de commande ./scriptIBD6.sh datas/brugada7chr22 datas/brugada8chr22.map datas/nvtsfdesirchr22.txt 
